#include "header.h"

/* -------------------------------------------------------------------------- */
void start(struct node* p){
/* -------------------------------------------------------------------------- */
double residue,res1;
int i,j,iter,ISXV_RUN=NO;
void find_dt(struct node*); 
void store_old(struct node*);
void flux_derv(int,struct node*,double**,double*);
void flux_derv_lim(int,struct node*,double**,double*);
void iterate(int,int,struct node*,void (*FLUX)(int,struct node*,double**,double*));
void find_residue(int,struct node*,double*,double*);
FILE* fres;

/* Coefficients of time-integration scheme */
AAA[0] = 0.0; AAA[1] = 3.0/4.0; AAA[2] = 1.0/3.0;
BBB[0] = 1.0; BBB[1] = 1.0/4.0; BBB[2] = 2.0/3.0;

fres = fopen(ResFile,"w");
residue  = 100.0;
iter     = LAST_ITER;
printf("Beginning of iterations\n");
time_elapsed = 0.0;

output(p);

var_old = (double **)calloc(no_grid_points,sizeof(double*));
if(var_old == NULL){ printf("var_old memory allocation failure.\n"); exit(0); }

if(time_mode == LOCAL){
	for(i=0;i<no_grid_points;i++)
    	var_old[i] = p[i].var;
}else {
	for(i=0;i<no_grid_points;i++)
		var_old[i] = (double *)calloc(5,sizeof(double));
}

for(i=0;i<no_grid_points;i++)
  for(j=0;j<5;j++)
    p[i].qx[j] = p[i].qy[j] = p[i].qz[j] = 0.0;


while(iter <= MAX_ITER && residue > MIN_RESIDUE){
	find_dt(p);
	store_old(p);
	for(i=0;i<rk_stage;i++) iterate(iter,i,p, 
        (void (*)(int,struct node*,double**,double*))(qmode ? flux_derv : flux_derv_lim));
	time_elapsed += p[0].dt;
	find_residue(iter,p,&residue,&res1);
	fprintf(fres,"%8d %15.6e\n",iter,residue);
	++iter;
	if(iter%ResFreq == 0){ 
		output(p); read_input(); parameters_used();
		}
	}

fclose(fres);

if(time_mode == GLOBAL)
	for(i=0;i<no_grid_points;i++) free(var_old[i]);
free(var_old);

}

/* -------------------------------------------------------------------------- */
void find_dt(struct node* p){ 
/* -------------------------------------------------------------------------- */
int i;
double DTIME;
double dt_node(int,struct node*);
double global_timestep(struct node*);

for(i=0;i<no_grid_points;i++)
	p[i].dt = cfl*dt_node(i,p);

if(time_mode == GLOBAL){ 
	DTIME = global_timestep(p);
	if(DTIME == 0.0){
		printf("Global time-step is zero.\n");
		exit(0);
		}
	for(i=0;i<no_grid_points;i++)
		p[i].dt = DTIME;
	}

}
/* -------------------------------------------------------------------------- */
double dt_node(int point,struct node* p){ /* Finds local time-step */
/* -------------------------------------------------------------------------- */
double dt,vel;
vel = sqrt( p[point].var[1]*p[point].var[1] + p[point].var[2]*p[point].var[2] 
            + p[point].var[3]*p[point].var[3]);
dt = p[point].dsm/(vel + 3.0*sqrt(p[point].var[4]/p[point].var[0]));
return dt;
}
/* -------------------------------------------------------------------------- */
double global_timestep(struct node* p){ /* Finds global time-step */
/* -------------------------------------------------------------------------- */
int i;
double DTIME;

DTIME = 1000.0;
for(i=0;i<no_grid_points;i++){
	if(p[i].dt < DTIME) DTIME = p[i].dt;
	}
return DTIME;
}
/* -------------------------------------------------------------------------- */
void store_old(struct node* p){
/* -------------------------------------------------------------------------- */
int i,j;

for(i=0;i<no_grid_points;i++)
	for(j=0;j<5;j++)
		p[i].var_old[j] = p[i].var[j];

}
/* -------------------------------------------------------------------------- */
void find_residue(int iter,struct node* p,double* res,double* res1){
/* -------------------------------------------------------------------------- */
int i,i_max_change;
double density_change;
double max_change,x_max_change,y_max_change, z_max_change;

max_change   = 0.0;
x_max_change = -1000;
y_max_change = -1000;
z_max_change = -1000;
i_max_change = -1;
*res         = 0.0;

for(i=0;i<no_grid_points;i++){
	density_change  = fabs(p[i].var[0] - p[i].var_old[0]);
	(*res)         += (density_change*density_change)/(p[i].var_old[0]*p[i].var_old[0]);
	if(density_change >= max_change){
		max_change   = density_change; 
		x_max_change = p[i].x;
		y_max_change = p[i].y;
		z_max_change = p[i].z;
		i_max_change = i;
		}
	}
*res = sqrt( (*res)/no_grid_points );
if(iter == LAST_ITER) (*res1) = (*res);
if( *res1 != 0.0) (*res) = (*res)/(*res1);
if(time_mode == GLOBAL)
printf("%5d %14.4e %14.4e %14.4e %10.5f %6d %3d %10.4f %10.4f %10.4f\n",iter,p[0].dt,(*res),
      time_elapsed,max_change,i_max_change,p[i_max_change].point_type,x_max_change,y_max_change,z_max_change);
else
printf("%6d %12.4e %10.5f %6d %3d %9.4f %7.4f %7.4f %7.4f %7.4f\n",iter,(*res),max_change,
	  i_max_change,p[i_max_change].point_type,x_max_change,y_max_change, z_max_change, Cl,Cd);
fflush(NULL);
}
